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Abstract 

In  this  paper,  a  dynamic  model  of  a  single  tubular  solid  oxide  fuel  cell  (SOFC)  unit  is  developed  using  the  control  volume  (CV)  approach.  The 
heat  transfer,  species  transportation,  and  electrochemical  reaction  effects  are  taken  into  account  in  a  collective  manner.  Using  this  model,  we 
study  the  spatial  distributions  of  a  series  of  state  variables  under  both  steady-state  and  transient  operations  and  evaluate  the  system  dynamic 
behavior.  The  analysis  shows  that  there  exists  non-uniform  current  contribution  and  Nernst  potential  distribution  along  the  longitudinal 
direction,  caused  by  the  non-uniform  fuel/gas  partial  pressures  along  the  flow  direction  and  the  temperature  distribution  in  the  anode/cathode 
channel  and  cell  body,  for  example.  The  classical  Nernst  potential  relation  is  revised  to  capture  the  characteristics  of  the  fuel  cell  under  time 
varying  external  load  voltage.  Comprehensive  numerical  simulations  are  carried  out  to  explore  the  underlying  dynamic  properties  in  typical 
SOFC  operations.  The  numerical  study  is  also  correlated  to  experimental  results  such  as  the  polarization  curve  and  power  density,  which 
shows  good  agreement. 
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1.  Introduction 

Solid  oxide  fuel  cells  (SOFCs)  are  very  promising  types 
of  high  temperature  fuel  cells  currently  being  considered 
for  commercial  power  source  applications.  Since  the  elec¬ 
trolyte  is  a  layer  of  ceramic  material  with  high  temperature- 
endurable  porous  media  electrodes,  the  SOFC  can  generally 
operate  at  high  temperature  range  (800-1000  °C).  High  oper¬ 
ating  temperature  can  lead  to  some  advantages,  such  as  high 
oxide-ion  conductivity  of  the  solid  oxide  electrolyte,  high 
rate  of  reaction  kinetics  and  hence  high-energy  conversion 
efficiency,  flexibility  of  fuel  type  due  to  reforming  capability 
of  hydrocarbons,  and  the  high  temperature  of  the  exiting  gas 
which  could  be  recycled  for  other  applications.  On  the  other 
hand,  the  high  temperature  operating  condition  also  imposes 


*  Corresponding  author.  Tel.:  +1  860  486  5911;  fax:  +1  860  486  5088. 
E-mail  address:  jtang@engr.uconn.edu  (J.  Tang). 

0378-7753/$  -  see  front  matter  ©  2004  Elsevier  B.V.  All  rights  reserved, 
doi:  10. 10 16/j.jpowsour.2004. 11.023 


some  challenges  which  include  potential  thermal  fatigue  fail¬ 
ure  of  the  cell  material,  sealing  under  high  temperature,  and 
the  associated  control  mechanism/power  required  to  maintain 
the  operation. 

Planar  and  tubular  geometries  are  the  two  most  popular 
forms  of  SOFC.  The  tubular  type  SOFC  is  relatively  ma¬ 
ture  in  terms  of  design  and  sealing,  while  the  planar  type 
shows  better  manufacturability.  The  core  of  a  SOFC  con¬ 
sists  of  a  porous  cathode  electrode,  electrolyte,  and  porous 
anode  electrode  assembly  (PEN),  where  the  electrochemical 
reaction  takes  place  and  electrical  energy  and  heat  are  gen¬ 
erated  directly  from  the  species  chemical  energy  [1,2].  The 
environmental  condition  of  the  electrochemical  reaction  has 
direct  effects  on  the  efficiency  of  the  fuel  cell.  Besides  the 
physical  properties  of  the  PEN,  the  cell  temperature,  partial 
pressures  of  hydrogen,  oxygen  and  vapor,  etc.  all  play  impor¬ 
tant  roles  in  the  electrochemical  reaction.  In  order  to  analyze 
the  complicated  interactions  between  the  various  factors  and 
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Fig.  1.  Tubular  cell  configuration. 


to  optimize  the  system  performance,  there  has  been  signif¬ 
icant  interest  in  developing  mathematical  models  for  SOFC 
[3-9].  Li  and  Chyu  [3]  proposed  a  2D  steady-state  model 
for  a  tubular  SOFC  in  a  stack  by  assuming  that  heat  and 
species  are  not  exchanged  between  a  cell  and  its  neighbors. 
The  model  did  not  consider  the  radiation  effect  among  cells. 
Similarly,  Nagata  et  al.  [4]  investigated  output  characteristics 
of  a  ID  steady-state  tubular  type  SOFC  whereas  the  radiation 
effect  is  included.  Haynes  [5]  developed  a  transient  model  to 
capture  the  heat  transfer  effects.  However,  some  conditions 
assumed  in  that  research  were  not  consistent  with  the  essen¬ 
tial  characteristics  of  the  transient  behavior.  For  example,  the 
inlet  gas/fuel  mole  fraction  was  assumed  constant.  Ota  et  al. 
[6]  explored  the  dynamic  behavior  of  a  micro-tube  SOFC. 
In  that  model,  they  imposed  the  incompressible  assumption 
of  gases  and  the  viscous  effect  was  not  considered.  Such  ef¬ 
fect,  nevertheless,  could  have  significant  impact  on  species 
partial  pressure  distribution  along  the  flow  direction  espe¬ 
cially  for  the  micro-tube,  and  in  turn  influence  the  voltage 
and  current  distributions  along  the  cell.  In  addition,  the  heat 
conduction  between  the  defined  slices  was  also  neglected  in 
that  model.  Cooper  et  al.  [7]  proposed  a  steady-state  model 
for  the  reaction  and  species  transport  process  of  a  tubular  type 
SOFC.  Since  only  steady-state  operation  is  considered,  the 
energy  equation  is  not  involved  in  the  modeling  and  analyti¬ 
cal  solution  with  asymptotic  expansion  can  be  obtained.  It  is 
worth  noting  that  the  electrochemical  characteristics  of  the 
cell  are  generally  treated  as  an  equivalent  circuit  model,  and 
the  Nernst  potential  is  usually  regarded  as  a  potential  source 
[3,4,6]. 

It  is  the  goal  of  this  study  to  develop  a  dynamic  model  of 
the  SOFC  that  is  capable  of  characterizing  the  transient/time- 
and  spatial-dependent  properties  of  critical  state  variables. 
Although  a  lot  of  work  has  been  dedicated  to  the  steady- state 
performance  evaluation  [3, 4, 7-9],  the  inclusion  of  transient 
dynamics  significantly  complicates  the  model  development. 
It  appears  that  the  direct  integration  of  time  dependency  into 
the  high  precision  steady-state  model  (such  as  the  one  given 
in  [3])  could  lead  to  insurmountable  computational  cost.  This 
is  especially  true  for  single  tubular  SOFCs  packed  in  a  stack. 
Clearly,  the  real  challenge  is  to  develop  a  model  that  can  be 
sufficiently  utilized  in  system  optimization  and  dynamic  con¬ 
trol  while  being  computationally  tractable.  In  this  research, 
building  upon  the  previous  studies  on  dynamic  modeling  of 
SOFC  and  using  the  control  volume  (CV)  approach,  we  de¬ 
velop  a  discrete  SOFC  model  which  incorporates  thoroughly 
the  effects  of  conduction,  convection,  radiation,  species  trans¬ 


portation,  and  electrochemical  reaction.  In  modeling  the  elec¬ 
trochemical  reaction,  we  use  a  revised  Nernst  potential  rela¬ 
tion  to  characterize  the  unique  flip-flop  behavior  in  a  SOFC 
under  time  varying  operating  conditions.  This  model  will  al¬ 
low  us  to  carry  out  detailed  analysis  on  current/voltage  dis¬ 
tributions  and  temperature  and,  for  example,  fuel/gas  partial 
pressure  distributions. 

2.  Modeling  of  SOFC 

2.7.  Single  tubular  SOFC  configuration 

We  consider  an  anode- supported  single  tubular  fuel  cell 
illustrated  in  Fig.  1  [1].  The  internal  cylinder  is  the  cell 
body,  covered  by  a  layer  of  thermal  insulator.  The  fuel  and 
oxygen  flow  in  opposite  directions,  which  could  make  the 
heat/temperature  distribution  along  the  cell  more  uniform. 
Three  basic  SOFC  designs  have  been  explored  in  the  open 
literature,  the  electrolyte-supported,  the  cathode-  and  anode- 
supported  designs.  Since  the  electrolyte-supported  SOFC  ex¬ 
hibits  very  high  electrolyte  resistance  under  moderate  tem¬ 
perature,  the  operation  must  be  at  very  high  temperature 
(^1000  °C).  Electrode- supported  designs  are  better  suited 
for  lower  temperature  operation.  It  has  been  argued  that  an 
anode- supported  SOFC  could  have  better  cell  output  power 
[1,2].  Without  loss  of  generality,  in  this  paper  we  focus  our 
attention  on  the  anode- supported  SOFC.  This  modeling  strat¬ 
egy  can  be  easily  extended  to  other  SOFC  designs. 

2.2.  Modeling  of  heat  transfer  effects 

A  SOFC  involves  complicated  heat  transfer,  species  trans¬ 
portation,  and  electrochemical  effects  which  are  highly  inter¬ 
active.  In  order  to  make  the  model  tractable  while  capturing 
the  fundamental  dynamic  behavior,  we  reduce  the  problem  to 
a  quasi-2D  analysis  based  on  the  changes  along  the  gas  flow 
direction  using  the  control  volume  (CV)  approach  [4-6,8].  In 
other  words,  while  the  flow  field  within  each  channel  is  con¬ 
sidered  as  ID  respectively,  the  cross-sectional  effects  (e.g., 
heat  transfer,  species  transportation,  etc.)  between  the  two 
concentric  channels  are  fully  taken  into  account.  As  shown 
in  Fig.  2,  the  cell  is  divided  into  serial  segments  by  the  planes 
perpendicular  to  fuel  flow  direction.  Each  segment  includes 
four  CVs,  i.e.,  anode  CV,  cell  CV,  cathode  CV,  and  thermal 
insulator  CV.  Within  the  control  volume,  we  assume  uni¬ 
form  physical  properties.  It  should  be  noted  that  the  physical 
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Fig.  2.  Mass  and  energy  balance  CV  illustration. 


properties  in  different  CVs  are  not  necessarily  the  same.  The 
electrochemical  reactions  take  place  in  the  cell  (PEN).  Due  to 
the  high  operating  temperature,  all  the  water  produced  by  the 
electrochemical  reactions  is  assumed  to  be  in  the  vapor  state. 
Here  we  assume  an  ideal  gas  mixture  within  the  cell,  which 
means  that  each  mixture  component  behaves  as  an  ideal  gas 
as  if  it  were  alone  [10].  For  the  control  volumes  of  the  PEN 
and  insulator,  only  the  energy  balance  equations  are  needed 
because  these  two  control  volumes  are  made  of  solid  mate¬ 
rials,  while  the  control  volumes  of  gas/fuel  must  satisfy  the 
balance  equations  of  mass/species,  momentum  and  energy. 

2.2.7.  Mass/species  conservation  equation 

First  we  consider  the  fuel/gas  CV.  The  mass/species  bal¬ 
ance  equation  for  this  CV  satisfies  [11]: 

Vi-f  +  (YpuA)i_l  -  ( YpuA)i+l  =  msurf  (1) 

where  V  is  the  control  volume,  p  the  species  density  within 
the  CV,  Y  the  species  mass  fraction,  A  the  cross-section  area, 
u  the  bulk  gas  flow  speed,  mSUrf  the  net  species  mass  flux 
through  the  connecting  surface  between  neighboring  CVs, 
and  the  subscript  i  designates  the  ith  CV. 

For  the  ith  anode  channel  CV,  the  entering  species  from 
the  (i  —  l)th  CV  and  the  leaving  species  to  the  (i  +  l)th  CV 
are  hydrogen  and  vapor.  The  diffusion  species  include  the 
hydrogen  for  the  electrochemical  reaction  (diffusing  into  the 
anode  electrode)  and  vapor  (generated  by  the  electrochemi¬ 
cal  reaction  and  diffusing  out  from  anode  electrode).  For  the 
cathode  channel  CV,  the  flow  stream  species  include  oxy¬ 
gen  and  nitrogen  from  air  or  pure  oxygen.  The  oxygen  is  the 


only  involved  species  diffusing  into  the  cathode  electrode 
under  electrochemical  reaction.  These  can  be  characterized 
by  plugging  relevant  species  parameters  into  Eq.  (1),  which 
will  yield  a  set  of  equations  that  describes  the  mass/species 
conservation  within  each  CV. 


2.2.2.  Momentum  equation 

Generally,  the  flow  field  also  satisfies  the  momentum  equa¬ 
tion  [11]: 


t/3  (E"PkYk) 

'  dt 


+ 


(fffupkYkuA) 


(F uPkYkuA ) 


longitudinal 


For  the  anode  channel  CV,  the  entering/leaving  momentum 
is  caused  by  main  flow  stream  hydrogen  and  water  vapor.  As 
compared  to  the  total  mass  of  main  flow  stream,  the  mass  di¬ 
rectly  involved  in  the  electrochemical  reaction  is  small,  and 
thus  the  momentum  effect  on  the  main  flow  stream  caused 
by  such  electrochemical  reaction  is  negligible.  The  cathode 
channel  that  involves  the  oxygen  and  nitrogen  flows  is  mod¬ 
eled  using  a  similar  assumption.  Again,  when  plugging  the 
respective  species  parameters  into  Eq.  (2),  we  can  get  a  set 
of  equations  that  characterize  the  momentum  balance  within 
the  CV. 


2.2.3.  Energy  balance  equation 

Consider  Fig.  2.  The  energy  balance  equation  for  the 
fuel/gas  CV  (corresponding  to  the  anode  channel)  can  be 
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expressed  as  follows  [10,11]: 


'  dt 


•  1  —  ( (9a  1 

v^mass,m'/_i  ^^mass,out^/+i 


(Gconv)/  (Gmass,diff^/ 


(3a) 


where  e  is  the  specific  internal  energy  within  the  CV, 
G  mass,  in  and  2mass,out the  energy  transfer  rates  along  the  flow 
stream,  Gconv  the  convective  heat  transfer  rate  between  the 
gas  CV  and  the  PEN/thermal  insulator  CV,  and  <2^ass  diff 
is  the  heat  convection  rate  between  the  gas  CV  and  the 
PEN/insulator  CV  due  to  mass  diffusion.  The  radiation  heat 
transfer  between  the  gas  CV  and  the  solid  CV  is  not  consid¬ 
ered. 

For  the  anode  channel,  the  flow  stream  energy  transfer 
rates  between  the  neighboring  CVs  can  be  represented  as, 
respectively: 


The  convection  heat  transfer  rate  between  the  gas  CV  and 
solid  CV  is  given  as  [12]: 


(GconvX-  =  hiAc{[(7 g)(.  -  (Ts\]  (3d) 

where  h  is  the  convective  heat  transfer  coefficient,  Ac  the 
circumferential  area  of  the  C  V,  Tg  the  gas  temperature  within 
the  CV,  and  Ts  is  the  temperature  of  the  solid  CV  adjacent  to 
the  gas  CV. 

The  heat  transfer  due  to  the  mass  diffusion  effect  includes 
the  gas  diffusion  entering  and  leaving  the  PEN  [12]: 

(GmassV*  =  E  V~(^(7g  -  Ts).  (3e) 


where  M  is  the  species  molecular  mass,  n  the  number  of 
electrons  per  mole  of  the  kth  reactant  species,  I  the  cur¬ 
rent  contribution  of  the  corresponding  CV  segment,  F  the 
Faraday’s  constant,  and  cp  is  the  specific  heat  of  species. 
For  the  anode  channel  CV,  hydrogen  and  vapor  are  in¬ 
volved  in  Eq.  (3e),  while  only  oxygen  is  involved  in 
the  cathode  electrode  diffusion  and  electrochemical  reac¬ 
tion. 

The  energy  balance  in  the  CVs  of  PEN  and  insulator  is 
given  as 

c| rj~i 

w,^  =  (U-i  +  E  (Gconv)/  T  E  (Q  mass)/ 

""^(Grad)/  —  (£out)/+l  H-  (Ggen)j  (4a) 

where  £jn  and  Eoui  are  the  heat  conduction  effects  due  to  the 
adjacent  solid  CVs,  Gconv  and  Gmass  are  given  in  Eqs.  (3d) 
and  (3e),  Grad  the  radiation  heat  transfer  rate  between  the 


cell  and  the  thermal  insulator,  and  Ggen  is  the  heat  genera¬ 
tion  rate  within  the  CV  due  to  electrochemical  reaction  and 
current  flow.  It  should  be  noted  that  since  the  PEN  and  insu¬ 
lator  CVs  are  composed  of  solid  materials,  only  the  energy 
balance  equation  is  needed.  The  heat  conduction  effects  can 
be  quantified  as  follows  [11]: 

7/_i  -  Ti 

(E[n)i_i  =  kA-—~ - -  (4b) 

Ax 

Ti  —  TV  i 

(Eout)/+i  =kA  1  1+1  (4c) 

Ax 


where  k  is  the  material  conductivity,  and  Ax  is  the  segment 
length. 

In  this  study,  we  take  into  account  the  radiation  effect 
between  the  cell  CV  and  its  adjacent  thermal  insulator  CV. 
Such  radiation  heat  transfer  can  be  approximated  by  using 
the  infinite  length  concentric  cylinder  theory  [12]: 


(Grad)/ 


1  /  ^PEN  +  (rpENMnsX(lMns) 


(4d) 


where  a  is  the  Stefan-Boltzmann  constant,  tpen  and  r\ns  are 
the  outside  and  inside  diameters  of  the  PEN  and  insulator, 
respectively,  and  s  is  the  surface  emissivity. 

The  electrochemical  reaction  heat  only  occurs  in  the  PEN 
layers  [4,14]: 


(Ggen)^  —  4Th2 


A/Th2  +  a  +  Rq  +  Rc)i 


(4e) 


where  ATTh2  is  the  low  heating  value  of  hydrogen,  Ra  is  the 
anode  overpotential  resistance,  Rq  is  the  cathode  overpoten¬ 
tial  resistance,  and  Rq  is  the  ohmic  resistance. 

The  above  derivation  completes  the  energy  balance  char¬ 
acterization  of  the  anode  channel,  the  cell  body,  and  the  ther¬ 
mal  insulator.  For  the  cathode  channel,  the  general  descrip¬ 
tion  is  very  similar  to  that  of  the  anode  channel.  In  addition 
to  the  species  difference,  however,  it  should  be  noted  that  for 
the  cathode  channel,  due  to  the  reverse  in  flow  direction,  the 
order  of  CVs  in  Eqs.  (1)— (3)  should  also  be  reversed,  i.e.,  i  +  1 
is  the  upstream  and  i—  1  is  the  down  stream.  Also,  for  the 
cathode  channel  C  V,  the  convection  involves  that  between  the 
PEN  and  the  cathode  channel  and  that  between  the  insulator 
and  the  cathode  channel,  while  the  convection  for  the  anode 
channel  only  involves  that  between  the  PEN  and  the  anode 
channel. 


2.3.  Modeling  of  electrochemical  reaction 

A  SOFC  is  an  electrochemical  device  that  converts  the 
chemical  energy  into  electrical  energy  directly.  At  the  cathode 
side,  the  oxygen  diffuses  through  the  electrode  and  reaches 
the  electrode/electrolyte  interface,  where  the  oxygen  is  elec- 
trochemically  reduced  into  oxygen  ions  by  consuming  the 
electrons  transported  through  the  external  circuit.  The  oxy¬ 
gen  ion  is  transported  through  the  electrolyte  to  the  anode 
side.  The  hydrogen  in  the  anode  channel  diffuses  through 
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the  electrode  and  reaches  the  electrode/electrolyte  interface, 
where  the  oxygen  ions  react  with  hydrogen  to  form  water 
and  electrons.  Electricity  is  generated  by  the  formed  elec¬ 
trons  at  the  anode  side.  The  oxidation  of  the  hydrogen  yields 
the  Nernst  potential  equation  [2,12,13]: 


—  A  G° 
IF 


+ 


RT 

2F 


In 


/JH2  ff  \ 

Pnio  j 


where  AG°  is  the  variation  of  standard  state  Gibbs’  free  en¬ 
ergy  of  oxidation  reaction  of  the  hydrogen,  R  the  universal 
gas  constant,  T  the  temperature  at  the  reaction  site,  Pu2  and 
Pr2o  are  the  ratios  of  the  hydrogen  and  water  vapor  par¬ 
tial  pressures  at  the  anode  electrode/electrolyte  interface  over 
the  standard  atmosphere  pressure  respectively,  and  Po2  is  the 
oxygen  counterpart  at  the  cathode  electrode/electrolyte  inter¬ 
face. 

The  operating  cell  potential  is  generally  lower  than  the 
Nernst  potential,  which  is  typically  caused  by  the  ohmic, 
concentration  and  activation  polarizations.  Based  upon  ex¬ 
perimental  characterizations,  the  anode  and  cathode  equiva¬ 
lent  resistances  are  quantified  as  [13,14]: 


4  F 
RT 


Ua\ 

rt) 


4  F 
~RT 


where  Ua  and  Uq  are  the  activation  energy  of  the  anode  and 
cathode,  respectively,  and  k\  and  kc  are  the  pre-exponential 
factors  that  are  obtained  experimentally. 

With  the  CV  approach  as  a  basis,  the  equivalent  fuel  cell 
circuit  [3,4,6,14]  can  be  formed  as  shown  in  Fig.  3.  Each  seg¬ 
ment  of  the  fuel  cell  constitutes  a  sub-fuel  cell.  It  is  assumed 
that  the  neighboring  sub-fuel  cells  have  no  direct  current  flow 
between  each  other  through  either  the  electrolyte  or  the  an¬ 
ode/cathode  electrodes,  since  the  internal  resistance  between 
them  is  much  larger. 

Traditionally,  the  Nernst  potential  of  a  sub-fuel  cell  is  sim¬ 
ply  treated  as  a  potential/voltage  source  [3,4,6,14],  and  the 
current  and  voltage  in  the  circuit  is  thus  calculated  through 


Fig.  4.  Illustration  of  modified  fuel  cell  Nernst  potential. 


Kirchhoff’s  law.  It  is  worth  noting  that  a  fuel  cell  can  be 
treated  as  a  voltage  source  only  if  the  Nernst  potential  is 
accumulated  high  enough  so  that  it  can  deliver  the  exter¬ 
nal  load  potential  requested.  If  the  external  load  potential  is 
higher  than  the  fuel  cell  Nernst  potential,  the  fuel  cell  will 
have  to  switch  the  role  from  potential  source  to  an  energy 
storage  element  [15].  This  can  be  illustrated  by  the  simula¬ 
tion  result  of  a  single  tubular  SOFC  as  given  in  Fig.  4  (de¬ 
tailed  discussions  will  be  given  in  Section  3).  As  will  be 
explained  later,  the  Nernst  potential  distribution  of  a  SOFC 
decays  along  the  longitudinal  direction.  Suppose  the  out¬ 
put  voltage  stays  at  a  constant  value,  for  example,  0.9  V,  as 
shown  in  Fig.  4.  The  Nernst  potentials  at  some  segments  are 
lower  than  the  external  load  voltage.  In  this  case,  instead 
of  the  Nernst  potential  providing  the  electrical  energy,  the 
external  voltage  or  the  Nernst  potentials  at  other  segments 
of  the  SOFC  will  provide  charge  to  the  sub-cells  until  the 
relevant  potentials  are  high  enough  to  satisfy  the  external 
voltage  requirement.  In  the  current  case,  the  downstream 
Nernst  potentials  will  finally  reach  the  level  of  the  exter¬ 
nal  load  voltage  as  shown  in  Fig.  4.  Depending  on  the  lev¬ 
els  of  sub-cell  Nernst  potentials  and  the  external  load  volt¬ 
age,  the  role  of  sub-cells  could  flip  back  and  forth.  Under 
these  considerations,  in  this  research  we  consider  the  sub¬ 
fuel  cell  as  a  combined  effect  of  the  Nernst  potential  and 
capacitor: 


-^Nernst  — 


E\E  =  - 
dV 

V\ —  = 
d  t 


A  G°  RT 
+  —  In 


P^P, 


1/2 


H2^02 

h20 


E  >  Vout 

E  <  Kut 

(8) 


where  Tout  is  the  terminal  output  voltage,  E  the  original 
Nernst  potential  produced  by  the  electrochemical  reaction, 
R  the  total  resistance  consisting  of  concentration,  activation 
and  ohmic  resistances,  C  the  capacitance  formed  by  sand- 
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wiched  electrode  and  electrolyte  of  sub-fuel  cell,  and  V  the 
sub-fuel  cell  potential  if  the  original  Nernst  potential  is  less 
than  the  output  voltage. 

The  Kirchhoff’s  voltage  and  current  laws  are  then  em¬ 
ployed  to  calculate  the  current  and  other  related  variables,  i.e., 

£out  =  (^Nernst)/  —  (Ra  +  Ra  +  Rc ),  /,  (9a) 


/out  =  7  (9b) 

i 

where  E0 ut  and  70Ut  are  the  load  voltage  and  the  load  current, 
respectively,  (^Nernst)/  and  7/  are  the  sub-fuel  cell  Nernst 
potential  and  the  associated  current.  The  three  polarization 
resistances  are  calculated  based  on  the  instant  conditions  of 
the  corresponding  sub-cell. 


Table  1 

Parameters  used  in  simulation 


Symbol 

Meaning 

Value 

Reference 

L 

Fuel  cell  length 

200  mm 

n 

Segment  number 

10 

Ax 

Segment  length 

20  mm 

ra 

Inside  radius  of  the  cell 

3.1  mm 

Outside  radius  of  the  cell 

3.9  mm 

A 

Inside  radius  of  the  insulator 

7.0  mm 

Al 

Outside  radius  of  the  insulator 

9.0  mm 

Mh2 

Molecule  weight  of  hydrogen 

2.016gmol_1 

[13] 

cpH2 

Specific  heat  of  hydrogen  (constant  pressure) 

14821  Jkg-1  K-1 

[10] 

Hydrogen  gas  constant 

4124.3  J  kg-1  K-1 

[13] 

cm2 

Specific  heat  of  hydrogen  (constant  volume) 

10697  Jkg-1  K_1 

[10] 

Ph2 

Hydrogen  density  at  900  K 

0.02723  kgm-3 

[10] 

/Xr2 

Dynamic  viscosity  of  hydrogen 

18.78  x  10_6kgm-1  s_1 

[10] 

AH 

Low  heating  value  of  hydrogen 

1.196  x  10s  J 

Well  known 

Molecule  weight  of  water  vapor 

18.02  g  mol-1 

[13] 

CpW20 

Specific  heat  of  water  vapor  (constant  pressure) 

2186  Jkg-1  K_1 

[10] 

Rei2o 

Water  vapor  gas  constant 

188.5  Jkg-1  K_1 

[13] 

Cv  H20 

Specific  heat  of  water  vapor  (constant  volume) 

1997.5  Jkg-1  K_1 

[10] 

Ph2o 

Water  vapor  density  at  850  K 

0.2579  kgm-3 

[10] 

Pn2o 

Dynamic  viscosity  of  water  vapor 

29.69  x  10~6  kg  m-1  s-1 

[10] 

mo2 

Molecule  weight  of  oxygen 

32  g  mol-1 

[13] 

cp02 

Specific  heat  of  oxygen  (constant  pressure) 

988.1  Jkg-1  K-1 

[10] 

r02 

Oxygen  gas  constant 

259.8  Jkg-1  K-1 

[13] 

Cv02 

Specific  heat  of  oxygen  (constant  volume) 

728.3  Jkg-1  K_1 

[10] 

Po2 

Oxygen  density  at  550  K 

0.7096  kgm-3 

[10] 

P02 

Dynamic  viscosity  of  oxygen 

31.97  x  10-6 kgm-1  s_1 

[10] 

ha 

Convective  coefficient  in  anode  channel 

2987  Wm-2  K-1 

a 

he 

Convective  coefficient  in  cathode  channel 

1322.8  Wm“2  K_1 

a 

hoc 

Convective  coefficient  insulator/surrounding 

10  Wm-2  K-1 

[10] 

Ppen 

Average  density  of  PEN 

6337.3  kgm“3 

[12] 

CpPEN 

Average  specific  heat  of  PEN 

594.3  J  kg" ‘K-1 

[12] 

kpEN 

Average  thermal  conductivity  of  PEN 

2.53  Wm-1  K_1 

[12] 

SPEN 

Average  emissivity  of  PEN 

0.33 

[10] 

Pins 

Density  of  insulator 

480  kg  m-3 

[10] 

Cplns 

Specific  heat  of  insulator 

10471kg-1  K-1 

[10] 

kins 

Thermal  conductivity  of  insulator 

0.059  Wm"1  K-1 

[10] 

£lns 

Thermal  emissivity  of  insulator 

0.09 

[10] 

<7 

Stefen-Boltzmann  constant 

5.669  x  10-8  W m-2  K-4 

[10] 

Ps 

Fuel/gas  source  pressure 

3.0  x  105  Pa 

b 

Ts 

Fuel/gas  source  temperature 

1073  K 

b 

Too 

Surrounding  temperature 

303  K 

b 

Ka_in 

Inlet  flow  coefficient  of  anode  channel 

8.7641  x  10-10kgPa-1  s_1 

c 

Ka_out 

Outlet  flow  coefficient  of  anode  channel 

4.3821  x  10-7  kg  Pa-1  s-1 

c 

Kc_in 

Inlet  flow  coefficient  of  cathode  channel 

3.1626  x  10-8  kg  Pa-1  s-1 

c 

^c_out 

Outlet  flow  coefficient  of  cathode  channel 

1.5813  x  10-5  kg  Pa-1  s_1 

c 

Ra 

Fuel  cell  ohmic  resistance 

0.0257  Q 

[10] 

Ea 

Anode  activation  energy 

llOkJ  mol-1 

[10] 

EC 

Cathode  activation  energy 

160kJ  mol-1 

[10] 

kA 

Anode  pre- exponential  factor 

2.13  x  108  Am-2 

[10] 

kc 

Cathode  pre-exponential  factor 

1.49  x  1010  Am-2 

[10] 

a  The  coefficients  are  calculated  by  assuming  a  fully  developed  laminar  flow  at  constant  wall  temperature  [6]. 
b  Parameters  are  assumed. 
c  Parameters  are  assumed  and  adjustable. 
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3.  Results  and  discussion 

The  above  section  presents  a  complete  set  of  mathematical 
relations  governing  the  transient  dynamic  behavior  of  SOFC. 
Using  these  equations,  we  develop  a  Matlab  code  to  simu¬ 
late  the  single  tubular  fuel  cell  operation.  All  physical  and 
geometrical  parameters  used  in  the  simulation  are  given  in 
Table  1 .  The  cell  is  divided  into  ten  segments,  where  the  two 
end  segments  do  not  involve  the  fuel/gas  diffusion  and  elec¬ 
trochemical  reaction.  The  mass  flow  rate  is  estimated  using 
the  nozzle  equation  [16].  The  model  is  highly  nonlinear  with 
distributed  parameters.  We  first  evaluate  the  steady-state  so¬ 
lution  using  an  iterative  procedure  and  the  dynamic  response 
is  then  calculated  using  stiff  solver  ode  15s  with  calculated 
initial  conditions. 

3.1.  Correlation  with  steady-state  experimental  results 

We  start  by  verifying  the  dynamic  model  using  a  correla¬ 
tion  of  the  available  experimental  data.  A  single  tubular  SOFC 
tested  is  a  nickel  and  yttria- stabilized  zirconia  (Ni-YSZ) 
cermet  anode-supported  tube  with  a  thin  YSZ  electrolyte 
coat  and  a  thin  lanthanum  strontium  magnate  (LSM)  cathode 
[1].  The  dimensions  and  the  relevant  physical  properties  of 
the  cell  are  given  in  Table  1 .  In  the  experiment,  the  fuel  cell 
temperature  is  kept  at  850  °C  by  using  an  electrical  furnace 
coupled  with  a  temperature  control  unit.  The  fuel  flowing 
through  the  anode  channel  is  pure  hydrogen,  whereas  the 
oxygen  is  obtained  directly  from  the  air  to  facilitate  the  elec¬ 
trochemical  reaction.  In  this  simplified  experiment,  the  cell 
does  not  have  the  thermal  insulator  layer.  We  specify  the  load 
current  and  correspondingly  measure  the  output  voltage.  This 
process  is  repeated  for  current  density  from  0  to  0.36  A  cm-2 
and  the  experimental  polarization  curve  and  power  density 
curve  are  obtained.  In  order  to  carry  out  the  correlation 
between  the  experiment  and  the  simulation,  we  insert  the 
experimental  parameters  and  conditions  into  the  dynamic 
model.  The  insulator  in  the  model  is  removed  to  fit  the  ex¬ 
perimental  condition,  and  the  cathode  channel  now  is  treated 
as  ambient.  As  the  temperature  of  the  cell  is  maintained  as 
constant,  the  energy  balance  equation  of  the  cell  body  is  not 
involved  in  the  simulation.  Under  these  assumptions,  we  per¬ 
form  steady-state  numerical  simulation,  and  the  polarization 
curve  and  the  power  density  curve  comparisons  are  given  in 
Figs.  5  and  6.  It  can  be  clearly  observed  that  the  simulation 
results  have  a  consistent  trend  when  compared  to  experimen¬ 
tal  data  and  the  model  can  predict  the  steady-state  behavior 
under  the  specific  experimental  conditions  reasonably  well. 

3.2.  Steady -state  performance  analysis 

A  more  detailed  analysis  is  carried  out  to  explore  the 
underlying  dynamic  properties  under  typical  SOFC  opera¬ 
tions.  The  following  analysis  is  performed  on  a  complete  cell 
unit  where  the  geometric  and  material  properties  are  listed  in 
Table  1 .  It  should  be  noted  that  the  key  parameters  are  those 


Fig.  5.  Experiment  and  simulation  comparison  of  polarization  curve. 


used  in  the  experimental  study  mentioned  above,  whereas  the 
cell  now  has  an  insulating  layer  and  operates  under  standard 
conditions  (without  a  temperature  control  unit).  The  species 
parameters  are  selected  based  on  previous  studies  [10,12]. 

To  evaluate  the  external  load  voltage  effect  on  the  perfor¬ 
mance  of  the  tubular  fuel  cell,  a  series  of  step  load  voltages 
(four  levels,  0.5,  0.7,  0.9,  and  1.1  V)  are  employed  as  the  in¬ 
put  to  the  system,  and  the  corresponding  steady-state  results 
are  illustrated  in  Figs.  7-14.  Fig.  7  shows  the  Nernst  potential 
distribution  along  the  longitudinal  direction  of  the  fuel  cell. 
When  the  external  load  voltage  is  0.5  V,  the  Nernst  potentials 
of  all  sub-fuel  cells  are  higher  than  the  external  load  volt¬ 
age,  and  thus  the  entire  tubular  fuel  cell  contributes  current 
to  the  external  circuit  based  on  the  Kirchhoff ’s  voltage  law.  In 
other  words,  under  such  a  scenario  the  Nernst  potentials  of  all 
sub-cells  obey  the  first  equation  shown  in  Eq.  (8).  This  result 
can  be  validated  by  the  current  distribution  given  in  Fig.  8, 
where  it  can  be  seen  that  all  sub-cells  contribute  (positive) 
current.  When  the  load  voltage  increases  to  0.7  V,  the  Nernst 


Fig.  6.  Experiment  and  simulation  comparison  of  output  power. 
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Fig.  7.  Fuel  cell  Nernst  potential  distribution  along  the  longitudinal  direc¬ 
tion. 

potentials  at  the  downstream  (in  the  direction  of  fuel  flow) 
segments  of  the  tubular  fuel  cell  become  very  close  to  the 
load  voltage.  Consequently,  while  the  upstream  segments  of 
the  fuel  cell  still  contributes  current,  the  current  contributions 
from  the  down  stream  segments  are  close  to  zero,  as  shown 
in  Fig.  8.  This  becomes  more  obvious  as  the  load  voltage 
further  increases  to  0.9  and  1.1  V.  In  these  cases,  while  the 
Nernst  potentials  of  several  upstream  segments  are  greater 
than  the  load  voltage,  the  Nernst  potentials  of  the  majority 
of  segments  (at  downstream  direction)  are  not  large  enough 
as  compared  to  the  load  voltage.  Correspondingly,  the  cur¬ 
rent  contributions  from  almost  the  entire  tubular  fuel  cell  are 
zero.  If  the  external  load  voltage  increases  further,  the  fuel 
cell  will  not  be  able  to  generate  current  but  instead  behave 
like  an  energy  storage  element,  i.e.,  a  capacitor.  This  can  be 
demonstrated  by  plotting  the  effective  length  of  the  tubular 
fuel  cell  (the  length  of  the  part  of  fuel  cell  that  contributes 


Fig.  8.  Current  contributions  along  the  flow  direction  with  different  load 
voltage. 


Load  Voltage  (V) 


Fig.  9.  Effective  length  of  fuel  cell  for  different  load  voltages. 

current  to  the  external  circuit)  versus  the  external  load  voltage 
(Fig.  9).  For  this  specific  SOFC  unit,  when  the  external  load 
voltage  is  lower  than  0.71  V,  all  segments  of  the  fuel  cell  will 
be  able  to  contribute  current  (and  thus  electrical  energy)  to  the 
external  circuit.  As  the  external  load  voltage  increases  (over 
0.71  V),  the  effective  length  of  the  fuel  cell  decreases  with 
only  a  portion  of  the  cell  contributing  current  and  useful  en¬ 
ergy.  The  effective  length  reaches  zero  when  the  external  load 
voltage  exceeds  1 . 1  V.  Such  property  is  very  important  in  the 
optimization  of  a  SOFC  under  specific  loading  conditions. 

As  shown  in  Eq.  (5),  the  water  vapor  partial  pressure  plays 
an  important  role  in  the  Nernst  potential.  In  Fig.  10  we  plot 
the  relation  between  the  external  load  voltage  and  the  vapor 
partial  pressure  distribution  along  the  longitudinal  direction. 
Generally,  the  vapor  partial  pressure  increases  when  the  exter¬ 
nal  load  voltage  decreases.  Meanwhile,  the  vapor  partial  pres¬ 
sure  increases  along  the  flow  direction.  Indeed,  lower  external 
load  voltage  corresponds  to  higher  current  contribution  (as 


Fig.  10.  Vapor  partial  pressure  distribution  along  the  flow  direction. 
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Cell  temperature  distribution 


Fig.  1 1 .  Cell  body  temperature  distribution  along  the  longitudinal  direction. 

illustrated  in  Fig.  8).  As  a  result,  the  hydrogen/oxygen  con¬ 
sumption  increases  and  thus  the  vapor  production  increases, 
yielding  higher  vapor  partial  pressure.  In  the  mean  time,  the 
heat  produced  by  the  electrochemical  reaction  and  the  three 
polarization  resistances  also  increase,  which  lead  to  an  in¬ 
crease  in  the  fuel  cell  body  temperature.  Fig.  11  shows  that 
the  fuel  cell  temperature  increases  with  a  decrease  of  the 
load  voltage.  The  increase  of  the  cell  body  temperature,  on 
the  other  hand,  leads  to  the  gas  temperature  increase  along 
the  flow  direction  for  both  anode  and  cathode  channels,  as 
shown  in  Figs.  12  and  13.  It  should  be  noted  that  the  change 
in  the  temperature  is  generally  insignificant.  One  interesting 


Fig.  13.  Cathode  channel  temperature  distribution  along  the  flow  direction. 

observation  is  that  the  change  of  gas  temperature  in  the  anode 
channel  is  more  significant  than  that  in  the  cathode  channel. 
The  reason  is  that  while  in  the  cathode  channel  both  oxygen 
diffusion  and  heat  convection  occur,  in  the  anode  channel, 
besides  the  hydrogen  diffusion  and  heat  convection,  the  dif¬ 
fusion  of  the  produced  vapor  also  takes  place.  In  the  case  that 
the  PEN  temperature  is  higher  than  that  of  the  gases  in  the 
channels,  the  vapor  diffusion  will  enforce  the  heat  transfer 
between  the  PEN  and  the  gases.  Fig.  14  shows  that  the  tem¬ 
perature  distribution  of  the  insulator  is  similar  to  those  of  the 
gases  in  the  cathode  channel.  To  illustrate  the  heat  transfer 
effects  to  the  PEN,  the  average  heat  flow  rates  towards  the 
PEN  for  load  cases  0.4  and  0.9  V  are  listed  in  Table  2.  The 


Fig.  12.  Anode  channel  temperature  distribution  along  the  flow  direction. 


Fig.  14.  Insulator  temperature  distribution  along  the  flow  direction. 


Table  2 


The  average  heat  transfer  to  the  PEN  (W) 


Generation 

Radiation 

Anode/convection 

Cathode/convection 

Anode/diffusion 

Cathode/diffusion 

0.4  V 

61.6625 

-2.6101 

-81.0218 

-326.1931 

-2.3537 

3.8833 

0.9  V 

0.0134 

-2.0952 

94.3834 

205.7795 

-0.0006 

0.0009 
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Fig.  15.  Voltage  contributions  transient  response. 


radiation  and  heat  transfer  due  to  mass  diffusion  are  negli¬ 
gible  as  compared  to  the  convective  heat  transfer.  When  the 
load  voltage  goes  up  from  0.4  to  0.9  V,  the  heat  generation 
due  to  electrochemical  reaction  and  the  equivalent  resistance 
decrease  because  the  current  contribution  goes  down,  which 
decrease  the  PEN  temperature  and  thus  the  heat  transfer  re¬ 
verses  the  direction. 

3.3.  Transient  behavior  analysis 

Finally  an  investigation  of  the  system  transient  behavior  is 
studied.  Such  transient  behaviors  are  critical  to  SOFC  control 
development.  In  general,  the  external  load  can  be  either  load 
current  or  load  voltage.  In  this  study,  we  use  external  step  volt¬ 
age  change  (step  load  voltage  experiencing  sudden  change 
from  0.6  to  0.5  V)  and  analyze  the  systems  overall  transient 
response.  The  inlet/outlet  gas/fuel  flow  rates  are  controlled 
by  using  the  traditional  nozzle  equation  [16].  The  flow  rate 


Fig.  16.  Current  contributions  transient  response. 


Fig.  17.  Inlet/outlet  hydrogen  flow  rate. 


coefficient  of  the  nozzle  usually  depends  on  the  temperature, 
pressure  difference  through  the  nozzle,  flow  rate  itself,  etc. 
To  simplify  the  simulation,  here  we  use  constant  flow  rate 
coefficients  (Table  1).  The  load  goes  down  from  0.6  to  0.5  V 
after  500s. 

In  Fig.  15,  E2-E9  curves  represent  the  time  history  of  the 
Nernst  potentials  of  the  segmented  sub-fuel  cells.  As  the  ex¬ 
ternal  load  voltage  goes  down  from  0.6  to  0.5  V,  the  Nernst 
potentials  of  all  sub-fuel  cells  go  down  and  show  a  slow  tran¬ 
sient  part  (about  1  h  in  the  current  simulation)  with  different 
magnitudes,  and  then  go  back  to  new  steady- state  values. 
The  transient  response  in  the  Nernst  potential  is  very  com¬ 
plex,  which  is  related  to  all  parametric  changes  associated 
with  the  entire  fuel  cell  activity.  Since  the  Nernst  potentials 
of  all  sub-cells  are  higher  than  the  external  load  voltages, 
all  sub-cells  contribute  current,  as  shown  in  Fig.  16.  Gener¬ 
ally,  these  results  are  in  total  agreement  with  the  steady- state 
observations. 


Fig.  18.  Vapor  partial  pressure  distribution  transient  response. 
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Fig.  19.  Cell  temperature  distribution  transient  response. 

After  the  load  voltage  experiences  sudden  change,  there 
are  overshoots  in  current  contributions  from  sub-cells. 
Moreover,  there  are  also  (delayed)  overshoots  in  both 
hydrogen  consumption  and  vapor  production,  as  shown 
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Fig.  20.  Anode  channel  temperature  distribution  transient  response. 
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Fig.  22.  Fuel  cell  output  power. 

in  Figs.  17  and  18.  The  increase  of  vapor  partial  pressure 
reduces  the  inlet  pressure  difference,  thus  the  inlet  hydrogen 
flow  rate  goes  down  for  a  period  of  time.  The  decrease  of 
inlet  hydrogen  flow  rate  and  the  increase  of  the  hydrogen 
consumption  will  reduce  the  outlet  hydrogen  flow  rate,  as 
illustrated  in  Fig.  18.  The  peak  (overshoot)  of  the  current  also 
increases  the  heat  production  and  thus  the  PEN  temperature 
increases  (based  on  Eq.  (4e)),  as  shown  in  Fig.  19.  The 
increase  of  PEN  temperature,  on  the  other  hand,  increases 
the  temperatures  of  the  anode  channel  gas  and  the  insulator, 
shown  in  Figs.  20  and  21 .  In  Fig.  22  we  plot  the  output  power 
history  of  the  fuel  cell  under  step  external  load  change. 


4.  Conclusions 

In  this  paper,  a  dynamic  model  of  a  single  tubular  SOFC 
combining  heat,  species  transportation  and  electrochemical 
reaction  effects  is  developed  to  study  its  steady-state  and 
transient  behavior.  The  polarization  performance  and  out¬ 
put  power  density  are  determined  by  a  variety  of  interacting 
fuel  cell  state  variables  such  as  fuel/gas  partial  pressure  dis¬ 
tribution,  heat  transfer  between  cell  body  and  gases.  Due  to 
the  non-uniform  distribution  of  the  Nernst  potential  along  the 
longitudinal  direction,  the  external  load  voltage  plays  an  im¬ 
portant  rule  in  the  cell  performance.  When  the  Nernst  poten¬ 
tials  of  some  segments  of  the  fuel  cell  are  lower  than  the  exter¬ 
nal  load  voltage,  these  segments  will  not  be  able  to  contribute 
current.  A  series  of  experimental  and  numerical  studies  are 
carried  out  to  verify  the  model  and  explore  the  underlying 
dynamic  behavior.  This  model  can  be  readily  used  in  system 
optimization  and  dynamic  control. 
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Fig.  21.  Insulator  temperature  distribution  transient  response. 
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